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ABSTRACT 

A  theoretical  analysis  suggests  that  the  path  centroid  variable  in  Feynman  path  in¬ 
tegration  occupies  a  central  role  in  the  behavior  of  the  real  time  position  autocorrelation 
function.  Based  on  this  analysis,  an  intriguing  quasiclassical  perspective  on  quantum  cor¬ 
relation  functions  emerges. 
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One  of  the  most  challenging  problems  in  condensed  matter  theory  is  the  computation 
of  quantum  time  correlation  functions.  The  satisfactory  solution  of  this  problem  will 
require  a  theoretical  formalism  which,  to  within  an  acceptable  degree  of  accuracy,  can 
be  implemented  efficiently  on  a  computer  for  general  many-body  systems.  The  present 
Commimication  is  devoted  to  one  such  formalism  which  is  based  on  the  properties  of  the 
path  centroid  variable  in  Feynman  path  integration^  which  is  defined  within  the  context 
of  the  imaginary  time  path  integral  as 

1 

In  his  derivation  of  an  approximate  quasiclassical  expression  for  the  partition  function,^ 
Feynman  recognized  that  the  path  centroid  was  a  particularly  useful  classical-like  variable 
about  which  to  define  a  quantmn  density,  often  called  the  '‘centroid  density” ,  defined  as 

pM  =  /■■■/  ^‘?(T)%c-qo)exp{-5[g(r)]/a}  ,  (2) 

where  5[9(t)]  is  the  imaginary  time  action  functional.  [The  quantum  partition  function 
can  be  obtained  from  Eq.  (2)  by  integration  over  the  position  of  the  centroid.  Qc-]  Feynman 
determined  a  variational  effective  potential  for  the  centroid  variable  using  the  approxima¬ 
tion  of  fi'ee  particle  motion  to  describe  the  action  for  the  imaginary  time  path  fluctuations 
about  the  centroid.  Since  Feynman’s  seminal  work,  several  authors^  have  improved  upon 
his  theory  by  employing  quadratic  variational  action  functionals  which  more  accurately 
approximate  the  exact  action  for  the  imaginary  time  path  fluctuations  about  the  centroid 
variable.  The  latter  approach,  therefore,  improves  the  accuracy  of  the  resulting  variational 
partition  function.  More  recently,  the  path  integral  centroid  perspective  has  been  extended 
to  calculate  equilibrium  averages  {A),  the  central  result  being  that  the  quantum  operator 
A  is  replaced  by  a  Gaussian  averaged  effective  classical  function.^  The  centroid  density 
has  also  become  the  centerpiece  of  a  path  integral  formulation  of  quantum  mechanical 
transition  state  theory.'*  To  our  knowledge,  however,  there  has  been  little  progress  in  the 
formulation  of  quantmn  dynamics  from  the  path  centroid  perspective. 


* 


We  have  recently  developed  a  theory^  for  quantum  imaginary  time  correlation  func¬ 
tions  in  terms  of  the  centroid  variable  and  its  corresponding  density  fxmction  [Eqs.  (i)  and 
(2)].  Quite  accurate  analytic  expressions  are  obtained  through  the  use  of  general  quadratic 
action  functionals,  a  properly  formulated  diagrammatic  perturbation  theory,  and  suitable 
renormalization  techniques.  In  principle,  once  a  good  approximation  to  the  exact  imagi¬ 
nary  time  correlation  function  is  obtained,  the  real  time  correlation  function  can  then  be 
determined  by  the  inverse  Wick  rotation  r  — ►  it  (i.e.,  by  analytic  continuation).  In  the 
course  of  this  research,  however,  an  even  more  promising  perspective  has  begrm  to  emerge 
on  the  role  of  the  centroid  variable  in  the  theory  of  quantum  time  correlation  functions. 
This  perspective,  which  is  the  subject  of  the  present  Communication,  offers  a  computa¬ 
tionally  powerful  algorithm  to  calculate  quantum  time  correlation  functions  in  many-body 
systems. 

To  begin  the  analysis,  it  proves  particularly  useful  to  define  the  centroid-constrained 
imaginary  time  correlation  function® 


r  (  n  \  -  /  •  •  •  /  ~  go)(<?('^)  ~  go)(g(0)  -  go)  ^^P{-S[q{T)/h]} 

J^  ■  J'Dq{r)6{qc-qo)exp{-S[q{r)]/h} 

This  correlation  function,  which  describes  the  correlations  of  the  imaginary  time  paths 
about  the  centroid  variable,  is  related  to  the  usual  correlation  function  C(r)  =  {q{T)q{Q)) 
by  the  relationship 


C(r)  =  (Cc(r,9c) +9c)pc  >  (4) 

where  the  subscript  “pc”  denotes  averaging  with  the  normalized  centroid  distribution 
PciQc)/  J  dqcpciqc)  [cf.  Eq.  (2)]. 

By  definition,  the  centroid  variable  plays  a  central  role  in  the  behavior  of  the  centroid- 
constrained  correlation  function  in  Eq.  (3).  A  more  subtle  issue,  however,  is  the  role  of 
the  centroid  variable  in  the  real  time  quantum  correlation  function.  In  a  formal  sense, 
this  information  can  be  extracted  from  the  exact  centroid-constrained  correlation  function 
C'c(ti9c)  through  the  analytic  continuation  r  — ►  it.  In  a  practical  sense,  however,  such 
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a  procedure,  whether  analytic  or  numerical,  is  not  tractable  unless  there  is  some  prior 
simplification  of  the  problem.  Fortunately,  such  a  simplification  can  be  achieved  through 
the  variational  determination  of  an  effective  quadratic  action  fimctional.  given  for  a  single 
degree-of-freedom  by^’® 

Seff[q{T)]  =  ^  dr  |yq(r)2  +  Vef/iqc)  +  ,  (5) 

where  Veff{qo)  is  the  variational  effective  centroid  potential  and  u>c  is  the  variationally 
determined  effective  harmonic  frequency,  both  quantities  being  evaluated  at  a  given  cen¬ 
troid  position.  As  in  other  work,  the  latter  parameter  is  determined  with  the  help  of 
the  properly  formulated^’^  Gibbs-Bogoliubov-Feynman  (GBF)  variational  principle®  for 
the  centroid  density  in  Eq.  (2).  The  optimized  frequency  is  given  by  the  solution  of  the 
transcendental  equation 


TTuj^  J dqV"(qc-^q)exp{-q^/2ac)  ,  (6a) 

where  the  effective  thermal  width  squared  in  the  optimized  reference  system  for  a  particular 
position  of  the  path  centroid  is  given  by 


etc 


fi^c/2 


77i/?u;|  itanh{hl3(jJc/‘2) 


(66) 


With  the  above  quadratic  reference  system  in  hand,  the  imaginary  time  correlation  function 
in  Eq.  (4)  can  be  determined  auialytically  for  the  reference  system,  giving® 


C(r) 


1  /■  j  ^  ^  h  cosh[uic{T  -  fiB/2)] 
Z  J  Pc(qc)  2^^  3inh{hl3Qj2) 


(7) 


where  Z  =  f  pc{qc)dqc  is  the  partition  function.  The  combined  and  second  aind  third  term 
in  Eq.  (7)  can  be  shown®  to  be  a  representation  in  the  effective  quadratic  reference  system 
of  the  simple  consttint  (9}^, where  (g)  is  the  thermal  average  of  the  variable  q.  FVom  now 
on,  those  two  terms  will  be  replaced  by  (g)^. 


Based  on  the  expression  in  Eq.  (7),  it  is  now  straightforward  to  analytically  continue 
the  above  function  to  obtain  the  real  time  correlation  function,  the  real  part  of  which  is 
given  by 

C(i)  =  ■  (*> 

At  this  level  of  the  theory,  the  real  time  correlation  function  is  seen  to  be  the  superposition 
of  centroid  correlation  fimctions  for  effective  harmonic  oscillators  defined  at  different  cen¬ 
troid  positions  Qc-  In  turn,  each  centroid  correlation  function  is  weighted  by  the  centroid 
density  for  that  value  of  Qc-  The  above  expression  for  the  correlation  function  is  there- 
foie  a  quantum  mechanical  generalization  of  the  classical  “instantaneous  normal  mode” 
perspective  for  condensed  matter  systems.^ 

Up  to  this  point,  the  arguments  have  remained  within  the  limits  of  rigorous  mathe¬ 
matical  analysis.  On  the  other  hand,  if  the  bounds  of  rigor  are  “loosened”  a  bit,  one  can 
arrive  at  an  even  more  compelling  result.  In  order  to  proceed  in  such  a  fashion,  it  proves 
beneficial  to  first  introduce  another  real  time  centroid  correlation  function,  given  by 

^  J  dqcpc{qc):^^^^cosiu;ct)  +  {q)^  ,  (9) 

which  is  the  exact  analog  of  the  classical  correlation  function  for  an  effective  harmonic 
oscillator  with  the  centroid  “frequency”  u>c-  Due  to  the  properties  of  the  Dirac  delta 
function,  the  above  correlation  function  is  related®  to  the  one  in  Eq.  (8)  through  their 
Fourier  transforms,  i.e., 

■ 

By  noting  that  the  first  term  on  the  right-hand-side  of  Eq.  (9)  describes  the  correlation  of 
fluctuations  about  the  mean  value  {qc)pc  —  (9)>  Eq.  (9)  can  then  be  rewritten  as 


C*{t)  =  ^  J  dQcPciqc){Qc{t)qc{0))^,  , 


where  {qc{t)Qc{0))we  is  a  quasidassicaL  position  correlation  function  calculated  for  dynamics 
on  the  locally  quadratic  potential  energy  surface  defined  for  each  initial  centroid  coordinate 


t 


•  • 


through  the  GBF  variational  procedure.  The  symbol  ■  ■)q^  denotes  initial  condition 
averaging  using  the  quadratic  approximation  to  the  centroid  density. 

At  this  point,  however,  one  might  ask  whether  Eq.  (11)  is.  in  fact,  an  approximation 
to  some  more  accurate  expression.  Some  insight  into  the  answer  is  provided  by  an  ex¬ 
amination  of  the  limiting  behavior  of  Eq.  (11)  for;  (a)  globally  harmonic  potentials,  and 
(b)  exact  classical  dynamics.  In  the  case  of  a  global  harmonic  potential,  the  correlation 
fimction  (oc(09c(0)}^^  in  Eq.  (11)  is  the  same  everywhere  [i.e.,  independent  of  the  in 
Eq.  (11)].  Since  f  dqcPciQc)  equals  unity,  Eq.  (11)  then  reduces  to  the  correlation 
function  (gc(^)9c(0))/),:-  In  the  global  harmonic  limit,  therefore,  one  can  simply  propagate 
the  centroid  variable  classically  using  forces  derived  from  the  effective  centroid  potential, 
average  the  initial  conditions  using  the  centroid  density,  and  calculate  the  classicaJ-like 
centroid  correlation  fimction 


c;{t)  =  (gc(09c(0)>p,  .  (12) 

The  exact  quantum  harmonic  correlation  function  is  then  determined  from  this  centroid 
correlation  function  via  Eq.  (10).  It  should  be  noted  that  for  general  nonlinear  potentials 
the  analytically  continued  effective  quadratic  result  in  Eq.  (9),  which  should  be  understood 
as  a  relatively  short  time  approximation  to  the  exact  correlation  function,  further  supports 
the  notion  of  a  quasiclassical  time  propagation  of  the  centroid. 

The  second  limiting  behavior  to  consider  is  the  case  of  exact  classical  dynamics.  The 
approximate  correlation  function  in  Eq.  (9)  is  given  classically  by 

Cciit)  =  J  dqd  PcliQci)  {qciit)qci{0))^,,  ,  (13) 

where  the  trajectories  qci{t)  are  generated  by  an  effective  quadratic  classical  potential 
with  frequency  uid-  The  distribution  based  on  this  potential  is  used  to  average  the  initial 
conditions  in  the  correlation  function,  while  the  second  level  of  averaging  is  performed  with 
the  exact  classical  density  Pdiqci)-  The  exact  classical  correlation  function  can  be  related  to 
Eq.  (13)  by  some  mental  mathematics.  In  particular,  if  the  potential  which  generates  the 
dynamics  and  averaging  in  the  correlation  function  {qd{t)<lci{0))Q^,  is  transformed  into  the 
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exact  potential,  that  correlation  function  becomes  the  exact  classical  correlation  function. 
The  "exterior  '  integration  over  classical  d^tribution  PciiQci)  then  becomes  inconsequential 
and  can  be  removed.  The  important  point  here  is  that  the  forces  which  generate  the 
classical  trajectories  in  (9d{f)9ci(0))  are  no  longer  the  approximate  linear  ones  as  in  Eq. 
(13),  but  are  instead  derived  from  the  same  nonlinear  potential  which  generates  the  exact 
equilibrium  distribution.  The  latter  distribution  is.  of  course,  the  classical  limit  of  the 
centroid  density.  One  might  then  argue  that  an  analogous  quantum  mechtmical  limit  of 
Eq.  (11)  is  one  in  which  the  forces  which  generate  the  centroid  trajectories  in  qdt)  are 
derived  from  the  same  effective  potential  which  gives  the  exact  centroid  distribution.  This 
notion  is.  of  course,  rigorously  correct  in  the  globally  harmonic  limit.  It  also  seems  logical 
that  the  forces  which  generate  the  centroid  trajectory  at  some  later  time  and  position  in 
space  should  be  no  different  from  the  forces  experienced  by  a  centroid  trajectory  which  is 
initiated  at  that  same  point  in  space. 

Based  on  the  above  arguments  and  Eqs.  (9)  —  (12),  the  following  picture  emerges  for 
a  centroid  density-based  calculation  of  the  position  time  correlation  function  in  quantum 
mechanics:  The  correlation  fimction  C*(t)  in  Eq.  (12)  should,  in  fact,  be  a  centroid 
correlation  fimction  which  is  calculated  using  the  exact  centroid  density  to  average  the 
initial  conditions.  The  centroid  trajectories  are  to  be  calculated  from  effective  classical 
equations  of  motion,  given  by 


mqc(t)  =  - 


dK(qc) 


where  the  effective  centroid  potential  is  given  by 

V'c(?c)  =  -fcBTlnlpciqc)]  ■  (15) 

The  centroid  force  Fdqc)  defined  in  the  right-hand-side  of  Eq.  (14)  is  a  kind  of  quantum 
mechanical  potentitil  of  mean  force,  given  explicitly  by 


Fdqc)  =  - 


/  ■  •  •  /  Vq{T)  6{qc  -  go)  ^[^(O)]  exp  {- S[q{T) / h]} 
/■••/  Vq{T)6{qc-qo)exp{-S[q{T)]/h} 
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To  determine  the  functional  form  of  the  centroid  momentum  distribution  used  to  average 
the  centroid  trajectory  initial  momenta,  one  needs  only  to  examine  the  action  functional 
for  the  imaginary  time  phase-space  path  integral,'*  given  by 


SE\p{T),q(T)]/h  = 


I.J  +  -  ip{T)g(T)| 

1 

hi  ■'’■''''’<’■>1 


Pn  ,  2n-T] 

2m  h0 


PnQ—n 


(17) 


where  {pn,Qn}  are  the  Fourier  modes  of  the  momentum  and  position  paths,  respectively. 
Remarkably,  the  momentum  centroid  pq  is  decoupled  from  position  coordinates  and  gen¬ 
erates  precisely  the  classical  momentum  distribution.  As  the  hnal  step  after  the  centroid 
correlation  function  C*{t)  in  Eq.  (12)  is  calculated,  the  real  part  of  the  quantum  position 
correlation  function  C{t)  can  be  determined  through  the  Fourier  transform  relationship  in 
Eq.  (10). 

The  above  "‘centroid  molecular  dynamics”  prescription  for  calculating  the  quantum 
position  correlation  function  has  both  rigorous  and  ad  hoc  elements  to  it,  though  it  is 
exact  for  globally  harmonic  systems  and  clearly  has  the  correct  classical  limit.  In  order 
to  subject  centroid  MD  to  a  numerical  test,  the  above  algorithm  was  applied  to  a  quartic 
oscillator  potential  given  by 


V{<l)  =  .  (18) 

with  the  par2imeters  g  =  2.4,  m  =  1.0,  h  =  1.0,  and  0  =  5.0.  The  exact  results  were 
obtained  by  finding  the  eigenvalues  and  eigenvectors  of  the  quartic  Hamiltonian  in  a  har¬ 
monic  oscillator  basis  set.  One  hundred  eigenstates  were  employed  to  yield  the  spectrum 
amd  time  correlation  function.  The  ground  state  energy  was  found  to  increase  by  12%  due 
to  the  quartic  perturbation  which  in  equivalent  to  a  large  anharmonic  shift  of  a  molec¬ 
ular  vibrational  spectnim.  The  centroid  MD  results  were  computed  using  the  centroid 
forces  from  the  variationally  optimized  quadratic  reference  system  in  Eq.  (5).  This  ref¬ 
erence  system  is  extremely  accurate  for  the  potential  in  Eq.  (18)  and  easy  to  implement 
numerically. 
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In  Fig.  1.  the  real  part  of  the  time  correlation  function  is  shown  for  the  exact  dynamics 
(stars),  the  centroid  MD  result  (solid  line),  the  analytically  continued  effective  harmonic 
result  in  Eq.  (8)  (dashed  line),  and  the  classical  limit  (dot-dashed  line).  The  centroid 
MD  is  clearly  the  superior  result,  capturing  the  shift  of  the  quantum  coherence  and  the 
(smaller)  dephasing  relative  to  the  classical  limit.  The  results  shown  in  Fig.  1  are  for  a 
very  low  temperature  system,  the  dynamics  of  which  are  dominated  by  only  the  ground 
and  first  excited  vibrational  states.  The  agreement  between  the  centroid  MD  and  the  exact 
results  will  become  even  better  at  higher  temperatmes. 

The  present  Communication  has  cont2uned  only  a  prefiminary  description  of  the  cen¬ 
troid  molecular  dynamics  method,  the  emphasis  being  on  the  position  time  correlation 
function.  Future  research  will  focus  on  an  extension  of  the  method  to  calculate  more  gen¬ 
eral  correlation  functions. a  more  rigorous  mathematical  justification  of  the  theory,^  the 
development  of  practical  algorithms  for  computing  centroid  dynamics,  and  apphcations 
of  centroid  MD  to  study  the  vibrational,  activated,  and  diffusional  dynamics  of  quantum 
particles  in  condensed  phase  system  .. 
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Figure  Caption 


Fig.  1;  A  plot  of  the  real  time  position  autocorrelation  function  for  the  quartic  oscillator 
described  in  Eq.  (18).  The  stars  show  the  numerically  exact  results,  the  sohd  line  is 
for  the  centroid  molecular  dynamics  method  (Eqs.  (12),  (14)-(17)|,  the  dashed  line  is 
for  the  analytically  continued  effective  qu2uiratic  theory  [Eq.  (8)],  and  the  dot-dashed 
line  is  for  classical  molecular  dynamics. 
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